!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine ucs (ro, mx, my, mz, bx, by, bz, en, nx, ng, flux_bnd, n_neighs, sys, dx, dt, gm)
!==================================================================================================================================|

   implicit none

   integer(4), intent(in) :: nx, ng, n_neighs, sys
   real(8), intent(in) :: dx, dt, gm
   real(8), dimension(n_neighs, sys), intent(inout) :: flux_bnd
   real(8), dimension(nx), intent(inout) :: ro, mx, my, mz, bx, by, bz, en

   integer(4) :: i, ip1, im1

   real(8), dimension(nx) :: ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre
   real(8), dimension(nx) :: ro_l, mx_l, my_l, mz_l, bx_l, by_l, bz_l, en_l
   real(8), dimension(nx) :: ro_r, mx_r, my_r, mz_r, bx_r, by_r, bz_r, en_r
   real(8), dimension(nx) :: ro_mod, mx_mod, my_mod, mz_mod, bx_mod, by_mod, bz_mod, en_mod
   real(8), dimension(nx) :: ro_g1, mx_g1, my_g1, mz_g1, bx_g1, by_g1, bz_g1, en_g1
   real(8), dimension(nx) :: ro_g2, mx_g2, my_g2, mz_g2, bx_g2, by_g2, bz_g2, en_g2
   real(8), dimension(nx) :: fro_pre, fmx_pre, fmy_pre, fmz_pre, fbx_pre, fby_pre, fbz_pre, fen_pre

   real(8) :: fro_l, fmx_l, fmy_l, fmz_l, fbx_l, fby_l, fbz_l, fen_l
   real(8) :: fro_r, fmx_r, fmy_r, fmz_r, fbx_r, fby_r, fbz_r, fen_r
   real(8) :: limiter
   integer(4) :: sgn
   real(8) :: dtdx, u1, u2, u3, vx, vy, vz, b2, v2, pr, gmm1, theta

!----------------------------------------------------------------------------------------------------------------------------------|
   dtdx = dt / dx
   gmm1 = gm - 1.0d0
   theta = 1.5

   do i = 2, nx - 1
      ip1 = i + 1
      im1 = i - 1

      u1 = (ro(i) - ro(im1)) * theta
      u2 = (ro(ip1) - ro(im1)) * 0.5d0
      u3 = (ro(ip1) - ro(i)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      ro_l(i) = ro(i) + 0.5d0 * limiter
      ro_r(i) = ro(i) - 0.5d0 * limiter
      ro_mod(i) = limiter

      u1 = (mx(i) - mx(im1)) * theta
      u2 = (mx(ip1) - mx(im1)) * 0.5d0
      u3 = (mx(ip1) - mx(i)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      mx_l(i) = mx(i) + 0.5d0 * limiter
      mx_r(i) = mx(i) - 0.5d0 * limiter
      mx_mod(i) = limiter

      u1 = (my(i) - my(im1)) * theta
      u2 = (my(ip1) - my(im1)) * 0.5d0
      u3 = (my(ip1) - my(i)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      my_l(i) = my(i) + 0.5d0 * limiter
      my_r(i) = my(i) - 0.5d0 * limiter
      my_mod(i) = limiter

      u1 = (mz(i) - mz(im1)) * theta
      u2 = (mz(ip1) - mz(im1)) * 0.5d0
      u3 = (mz(ip1) - mz(i)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      mz_l(i) = mz(i) + 0.5d0 * limiter
      mz_r(i) = mz(i) - 0.5d0 * limiter
      mz_mod(i) = limiter

      u1 = (bx(i) - bx(im1)) * theta
      u2 = (bx(ip1) - bx(im1)) * 0.5d0
      u3 = (bx(ip1) - bx(i)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      bx_l(i) = bx(i) + 0.5d0 * limiter
      bx_r(i) = bx(i) - 0.5d0 * limiter
      bx_mod(i) = limiter

      u1 = (by(i) - by(im1)) * theta
      u2 = (by(ip1) - by(im1)) * 0.5d0
      u3 = (by(ip1) - by(i)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      by_l(i) = by(i) + 0.5d0 * limiter
      by_r(i) = by(i) - 0.5d0 * limiter
      by_mod(i) = limiter

      u1 = (bz(i) - bz(im1)) * theta
      u2 = (bz(ip1) - bz(im1)) * 0.5d0
      u3 = (bz(ip1) - bz(i)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      bz_l(i) = bz(i) + 0.5d0 * limiter
      bz_r(i) = bz(i) - 0.5d0 * limiter
      bz_mod(i) = limiter

      u1 = (en(i) - en(im1)) * theta
      u2 = (en(ip1) - en(im1)) * 0.5d0
      u3 = (en(ip1) - en(i)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      en_l(i) = en(i) + 0.5d0 * limiter
      en_r(i) = en(i) - 0.5d0 * limiter
      en_mod(i) = limiter

   enddo

   do i = 1, nx
      vx = mx_l(i) / ro_l(i)
      vy = my_l(i) / ro_l(i)
      vz = mz_l(i) / ro_l(i)
      b2 = bx_l(i) * bx_l(i) + by_l(i) * by_l(i) + bz_l(i) * bz_l(i)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en_l(i) - v2 * ro_l(i) / 2.0d0 - b2 / 2.0d0) * gmm1
      fro_l = mx_l(i)
	   fmx_l = vx * mx_l(i) + pr + b2 / 2.0d0 - bx_l(i) * bx_l(i)
	   fmy_l = vx * my_l(i) - bx_l(i) * by_l(i)
	   fmz_l = vx * mz_l(i) - bx_l(i) * bz_l(i)
	   fbx_l = 0.0d0
	   fby_l = vx * by_l(i) - bx_l(i) * vy
	   fbz_l = vx * bz_l(i) - bx_l(i) * vz
	   fen_l = (en_l(i) + pr + b2 / 2.0d0) * vx - (bx_l(i) * vx + by_l(i) * vy + bz_l(i) * vz) * bx_l(i)

      vx = mx_r(i) / ro_r(i)
      vy = my_r(i) / ro_r(i)
      vz = mz_r(i) / ro_r(i)
      b2 = bx_r(i) * bx_r(i) + by_r(i) * by_r(i) + bz_r(i) * bz_r(i)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en_r(i) - v2 * ro_r(i) / 2.0d0 - b2 / 2.0d0) * gmm1
      fro_r = mx_r(i)
	   fmx_r = vx * mx_r(i) + pr + b2 / 2.0d0 - bx_r(i) * bx_r(i)
	   fmy_r = vx * my_r(i) - bx_r(i) * by_r(i)
	   fmz_r = vx * mz_r(i) - bx_r(i) * bz_r(i)
	   fbx_r = 0.0d0
	   fby_r = vx * by_r(i) - bx_r(i) * vy
	   fbz_r = vx * bz_r(i) - bx_r(i) * vz
	   fen_r = (en_r(i) + pr + b2 / 2.0d0) * vx - (bx_r(i) * vx + by_r(i) * vy + bz_r(i) * vz) * bx_r(i)

      ro_pre(i) = ro(i) - 0.5d0 * dtdx * (fro_l - fro_r)
      mx_pre(i) = mx(i) - 0.5d0 * dtdx * (fmx_l - fmx_r)
      my_pre(i) = my(i) - 0.5d0 * dtdx * (fmy_l - fmy_r)
      mz_pre(i) = mz(i) - 0.5d0 * dtdx * (fmz_l - fmz_r)
      bx_pre(i) = bx(i) - 0.5d0 * dtdx * (fbx_l - fbx_r)
      by_pre(i) = by(i) - 0.5d0 * dtdx * (fby_l - fby_r)
      bz_pre(i) = bz(i) - 0.5d0 * dtdx * (fbz_l - fbz_r)
      en_pre(i) = en(i) - 0.5d0 * dtdx * (fen_l - fen_r)

      vx = mx_pre(i) / ro_pre(i)
      vy = my_pre(i) / ro_pre(i)
      vz = mz_pre(i) / ro_pre(i)
      b2 = bx_pre(i) * bx_pre(i) + by_pre(i) * by_pre(i) + bz_pre(i) * bz_pre(i)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en_pre(i) - v2 * ro_pre(i) / 2.0d0 - b2 / 2.0d0) * gmm1
      fro_pre(i) = mx_pre(i)
	   fmx_pre(i) = vx * mx_pre(i) + pr + b2 / 2.0d0 - bx_pre(i) * bx_pre(i)
	   fmy_pre(i) = vx * my_pre(i) - bx_pre(i) * by_pre(i)
	   fmz_pre(i) = vx * mz_pre(i) - bx_pre(i) * bz_pre(i)
	   fbx_pre(i) = 0.0d0
	   fby_pre(i) = vx * by_pre(i) - bx_pre(i) * vy
	   fbz_pre(i) = vx * bz_pre(i) - bx_pre(i) * vz
	   fen_pre(i) = (en_pre(i) + pr + b2 / 2.0d0) * vx - (bx_pre(i) * vx + by_pre(i) * vy + bz_pre(i) * vz) * bx_pre(i)
   enddo

   do i = 1, nx - 1
      ip1 = i + 1
      ro_g1(i) = 0.5d0 * (ro(i) + ro(ip1)) + 0.125d0 * (ro_mod(i) - ro_mod(ip1)) - dtdx * (fro_pre(ip1) - fro_pre(i))
      mx_g1(i) = 0.5d0 * (mx(i) + mx(ip1)) + 0.125d0 * (mx_mod(i) - mx_mod(ip1)) - dtdx * (fmx_pre(ip1) - fmx_pre(i))
      my_g1(i) = 0.5d0 * (my(i) + my(ip1)) + 0.125d0 * (my_mod(i) - my_mod(ip1)) - dtdx * (fmy_pre(ip1) - fmy_pre(i))
      mz_g1(i) = 0.5d0 * (mz(i) + mz(ip1)) + 0.125d0 * (mz_mod(i) - mz_mod(ip1)) - dtdx * (fmz_pre(ip1) - fmz_pre(i))
      bx_g1(i) = 0.5d0 * (bx(i) + bx(ip1)) + 0.125d0 * (bx_mod(i) - bx_mod(ip1)) - dtdx * (fbx_pre(ip1) - fbx_pre(i))
      by_g1(i) = 0.5d0 * (by(i) + by(ip1)) + 0.125d0 * (by_mod(i) - by_mod(ip1)) - dtdx * (fby_pre(ip1) - fby_pre(i))
      bz_g1(i) = 0.5d0 * (bz(i) + bz(ip1)) + 0.125d0 * (bz_mod(i) - bz_mod(ip1)) - dtdx * (fbz_pre(ip1) - fbz_pre(i))
      en_g1(i) = 0.5d0 * (en(i) + en(ip1)) + 0.125d0 * (en_mod(i) - en_mod(ip1)) - dtdx * (fen_pre(ip1) - fen_pre(i))
   enddo

   do i = 2, nx - 1
      ip1 = i + 1
      im1 = i - 1

      u1 = (ro_g1(i) - ro_g1(im1)) * theta
      u2 = (ro_g1(ip1) - ro_g1(im1)) * 0.5d0
      u3 = (ro_g1(ip1) - ro_g1(i)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      ro_g2(i) = limiter

      u1 = (mx_g1(i) - mx_g1(im1)) * theta
      u2 = (mx_g1(ip1) - mx_g1(im1)) * 0.5d0
      u3 = (mx_g1(ip1) - mx_g1(i)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      mx_g2(i) = limiter

      u1 = (my_g1(i) - my_g1(im1)) * theta
      u2 = (my_g1(ip1) - my_g1(im1)) * 0.5d0
      u3 = (my_g1(ip1) - my_g1(i)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      my_g2(i) = limiter

      u1 = (mz_g1(i) - mz_g1(im1)) * theta
      u2 = (mz_g1(ip1) - mz_g1(im1)) * 0.5d0
      u3 = (mz_g1(ip1) - mz_g1(i)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      mz_g2(i) = limiter

      u1 = (bx_g1(i) - bx_g1(im1)) * theta
      u2 = (bx_g1(ip1) - bx_g1(im1)) * 0.5d0
      u3 = (bx_g1(ip1) - bx_g1(i)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      bx_g2(i) = limiter

      u1 = (by_g1(i) - by_g1(im1)) * theta
      u2 = (by_g1(ip1) - by_g1(im1)) * 0.5d0
      u3 = (by_g1(ip1) - by_g1(i)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      by_g2(i) = limiter

      u1 = (bz_g1(i) - bz_g1(im1)) * theta
      u2 = (bz_g1(ip1) - bz_g1(im1)) * 0.5d0
      u3 = (bz_g1(ip1) - bz_g1(i)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      bz_g2(i) = limiter

      u1 = (en_g1(i) - en_g1(im1)) * theta
      u2 = (en_g1(ip1) - en_g1(im1)) * 0.5d0
      u3 = (en_g1(ip1) - en_g1(i)) * theta
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (min (abs (u1), sgn * u2), sgn * u3))
      en_g2(i) = limiter
   enddo

   do i = 2, nx
      im1 = i - 1
      ro(i) = 0.5d0 * (ro_g1(im1) + ro_g1(i)) + 0.125d0 * (ro_g2(im1) - ro_g2(i))
      mx(i) = 0.5d0 * (mx_g1(im1) + mx_g1(i)) + 0.125d0 * (mx_g2(im1) - mx_g2(i))
      my(i) = 0.5d0 * (my_g1(im1) + my_g1(i)) + 0.125d0 * (my_g2(im1) - my_g2(i))
      mz(i) = 0.5d0 * (mz_g1(im1) + mz_g1(i)) + 0.125d0 * (mz_g2(im1) - mz_g2(i))
      bx(i) = 0.5d0 * (bx_g1(im1) + bx_g1(i)) + 0.125d0 * (bx_g2(im1) - bx_g2(i))
      by(i) = 0.5d0 * (by_g1(im1) + by_g1(i)) + 0.125d0 * (by_g2(im1) - by_g2(i))
      bz(i) = 0.5d0 * (bz_g1(im1) + bz_g1(i)) + 0.125d0 * (bz_g2(im1) - bz_g2(i))
      en(i) = 0.5d0 * (en_g1(im1) + en_g1(i)) + 0.125d0 * (en_g2(im1) - en_g2(i))
   enddo

   flux_bnd(1, 1) = fro_pre(ng)
   flux_bnd(2, 1) = fro_pre(nx - ng)
   flux_bnd(1, 2) = fmx_pre(ng)
   flux_bnd(2, 2) = fmx_pre(nx - ng)
   flux_bnd(1, 3) = fmy_pre(ng)
   flux_bnd(2, 3) = fmy_pre(nx - ng)
   flux_bnd(1, 4) = fmz_pre(ng)
   flux_bnd(2, 4) = fmz_pre(nx - ng)
   flux_bnd(1, 5) = fbx_pre(ng)
   flux_bnd(2, 5) = fbx_pre(nx - ng)
   flux_bnd(1, 6) = fby_pre(ng)
   flux_bnd(2, 6) = fby_pre(nx - ng)
   flux_bnd(1, 7) = fbz_pre(ng)
   flux_bnd(2, 7) = fbz_pre(nx - ng)
   flux_bnd(1, 8) = fen_pre(ng)
   flux_bnd(2, 8) = fen_pre(nx - ng)

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine ucs
